# This script makes all the plots based on GSS data (Figures 1, 2, 3, and C1 in the Appendix)

# Load packages and install them if necessary (automatic)
need = c(
  "dplyr",           #data structuring/cleaning
  "extrafont",       #plotting (fonts)
  "ggplot2",         #plotting
  "ggpubr",          #plotting (combining plots)
  "jtools",          #plotting (theme)
  "tidyr"            #data structuring/cleaning
)     
have = need %in% rownames(installed.packages())
if(any(!have)) {install.packages(need[!have])}
pack <- lapply(need, library, character.only = TRUE)
# Clean environment
rm(have, need, pack)

# Load data
load(file = "D:/dv_rep/Datasets/anes_boot_Plot.Rdata")

# Load fonts
extrafont::loadfonts()

#### FIG D1 (COHORT) ####
cohortDiv <- anes_Plot %>%
  filter(
    group == "Cohort" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  ylim(c(0.1, 0.30)) + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

cohortAlign <- anes_Plot %>%
  filter(
    group == "Cohort" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  scale_y_continuous(limits = c(0, 0.75), breaks = seq(0, 0.70, by = 0.10), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

legend_genMain <- ggpubr::get_legend(cohortDiv)

anes_gen_Main <- ggarrange(
  cohortDiv + theme(legend.position = "none"),
  cohortAlign + theme(legend.position = "none"),
  ncol = 2
)

anes_gen_Main <- ggarrange(
  anes_gen_Main,
  legend_genMain,
  ncol = 1,
  heights = c(20, 1)
)

anes_Main_genPlot <- annotate_figure(anes_gen_Main, top = text_grob("Cohort effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))

# Save plot as .png and as .pdf
ggsave(anes_Main_genPlot, file = "figd1.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

ggsave(anes_Main_genPlot, file = "figd1.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))


#### FIG D2 (PERIOD) ####

yearDiv <- anes_Plot %>%
  filter(
    group == "Year" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term)) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls, linetype = controls)) + 
  ylim(c(0.1, 0.30)) + 
  geom_point(position = position_dodge(width = 0.5), size = 2.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  #scale_x_continuous(breaks = seq(1950, 2021, by = 10), labels = seq(1950, 2021, by = 10)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), size = 3, linetype = c(1, 2))),
    shape = "none",
    linetype = "none"
  )

yearAlign <- anes_Plot %>%
  filter(
    group == "Year" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term)) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls, linetype = controls)) + 
  scale_y_continuous(limits = c(0, 0.85), breaks = seq(0, 0.80, by = 0.20), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 2.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), linetype = c(1, 2))),
    shape = "none",
    linetype = "none"
  )


legend_yearMain <- get_legend(yearDiv)

anes_year_Main <- ggarrange(
  yearDiv + theme(legend.position = "none"),
  yearAlign + theme(legend.position = "none"),
  ncol = 2
)

anes_year_Main <- ggarrange(
  anes_year_Main,
  legend_yearMain,
  ncol = 1,
  heights = c(20, 1)
)

anes_Main_yearPlot <- annotate_figure(anes_year_Main, top = text_grob("Period effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))

# Save output
ggsave(anes_Main_yearPlot, file = "figd2.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

ggsave(anes_Main_yearPlot, file = "figd2.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))


#### FIG D3 (AGE) ####

ageDiv <- anes_Plot %>%
  filter(
    group == "Age" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  ylim(c(0.15, 0.30)) + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), size = 3)),
    shape = "none",
    linetype = "none"
  )

ageAlign <- anes_Plot %>%
  filter(
    group == "Age" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  scale_y_continuous(limits = c(0, 0.70), breaks = seq(0, 0.70, by = 0.20), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

legend_ageMain <- get_legend(ageDiv)

anes_age_Main <- ggarrange(
  ageDiv + theme(legend.position = "none"),
  ageAlign + theme(legend.position = "none"),
  ncol = 2
)

anes_age_Main <- ggarrange(
  anes_age_Main,
  legend_ageMain,
  ncol = 1,
  heights = c(20, 1)
)

anes_Main_agePlot <- annotate_figure(anes_age_Main, top = text_grob("Age effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))

# Save output
ggsave(anes_Main_agePlot, file = "figd3.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

ggsave(anes_Main_agePlot, file = "figd3.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))


#### FIG C2 (MEAN LEVELS) ####

cohortMean <- anes_Plot %>%
  filter(group == "Cohort") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  mutate(term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Cohort") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

yearMean <- anes_Plot %>%
  filter(group == "Year") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Period") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

ageMean <- anes_Plot %>%
  filter(group == "Age") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  mutate(term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Age") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )


legend_Mean <- get_legend(yearMean)

mean_Plot_anes <- ggarrange(
  cohortMean + theme(legend.position = "none"),
  yearMean   + theme(legend.position = "none"),
  ageMean    + theme(legend.position = "none"),
  ncol = 3
)

mean_Plot_anes <- ggarrange(
  mean_Plot_anes,
  legend_Mean,
  ncol = 1,
  heights = c(20, 1)
)

mean_Plot_anes <- annotate_figure(mean_Plot_anes, top = text_grob("Mean-levels (ANES)", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))

# Save output
ggsave(mean_Plot_anes, file = "figc2.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

ggsave(mean_Plot_anes, file = "figc2.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))